clear
cd ..\..

local fname analyze_1_tables
capture log close
local c_date = c(current_date)
local c_time = c(current_time)
local c_time_date = "`c_date'"+"_"+"`c_time'"
local c_time_date = subinstr("`c_time_date '", ":", "_", .)
local c_time_date = subinstr("`c_time_date '", " ", "_", .)
log using "output/logs/`fname'_`c_time_date'.txt", replace text

********************************************************************************
*                                Define Macros                                 *
******************************************************************************** 
macro drop _all

* specify baseline consumption variable
global pre_avg "pre_avg_12"

* common covariates
global b $pre_avg // inherit baseline consumption variable from pre_avg
global c "income 	member 	uni_size 	building_size 	vintage	married"
global waveXmonth "month#opower_paper_wave"

* heterogeneous treatment effects - use relevant subset of demographic variables
global c_i "		member 	uni_size 	building_size 	vintage married"
global c_m "income 			uni_size 	building_size 	vintage married"
global c_u "income 	member 				building_size 	vintage married"
global c_b "income 	member 	uni_size 					vintage married"
global c_v "income 	member 	uni_size 	building_size 			married"

********************************************************************************
*                           Opower accounts by wave                            *
******************************************************************************** 
use .\intermediate_data\elec_2014-2018_mainfile, clear

keep ky_ba opower opower_paper_wave
duplicates drop
cap drop opower_unique
gen opower_unique = 1
cap drop opower_part_unique
gen opower_part_unique = 1 if opower==1
collapse (sum) opower_unique opower_part_unique, by (opower_paper_wave)

* add wave 5 and remove missing values
tsset opower_paper_wave
tsfill
replace opower_unique = 0 if opower_paper_wave==5
replace opower_part_unique = 0 if opower_paper_wave==5
drop if opower_paper_wave==.

estpost tabstat opower_unique opower_part_unique, by(opower_paper_wave) nototal
esttab . using ./output/tables/T2.tex, replace booktabs ///
	cells("opower_unique(fmt(%9.0fc)) opower_part_unique(fmt(%9.0fc))") ///
	nonumber noobs nomtitle ///
	collabels("Electric accounts in RCT" "Treated accounts", lhs("Wave"))

********************************************************************************
*                            Summary statistics                                *
******************************************************************************** 
use .\intermediate_data\elec_wave367_winterpre.dta, clear

eststo clear
quietly reg tot_kwh opower##post_opower $b $c // restrict to estimation sample

eststo: qui estpost sum tot_kwh income member building_size uni_size vintage married if e(sample), detail // add median to table

distinct ky_ba
estadd scalar hh = `r(ndistinct)' // get number of distinct households

esttab . using ./output/tables/T1.tex, replace booktabs ///
	cells("mean(fmt(%09.3gc)) p50(fmt(%09.3gc)) sd(fmt(%9.3gc)) min max") ///
	stats(hh N, fmt(%9.0fc) labels("Households" "Observations")) ///
	label nomtitle nonumber

********************************************************************************
*                            Opower treatment effect                           *
********************************************************************************
use .\intermediate_data\elec_wave367_winterpre.dta, clear

// all waves combined 
eststo clear
eststo: qui reghdfe tot_kwh opower##post_opower $b $c, absorb($waveXmonth) vce(cluster ky_ba)
qui sum tot_kwh if opower==0  // get pre-treatment mean for pooled sample
local r_mean_hat: di %9.0fc `r(mean)' // round the number 
estadd local control_mean `r_mean_hat' // store in local for display in table 

// wave 3, wave 6, wave 7 separate
foreach i of numlist 3 6 7 {
	eststo: qui reghdfe tot_kwh opower##post_opower $b $c if opower_paper_wave==`i', absorb($waveXmonth) vce(cluster ky_ba)
	qui sum tot_kwh if opower==0 & opower_paper_wave==`i'
	local r_mean_hat: di %9.0fc `r(mean)'
	estadd local control_mean `r_mean_hat'
}

// add FE indicators 
estfe . est*, labels($waveXmonth "Wave $\times$ year-month FE")
return list // create local used in table creation

// Export to tex
esttab using ./output/tables/T3.tex, replace booktabs ///
	b(%9.2f) se(%9.2f) ///
	star(* 0.1 ** 0.05 *** 0.01)  ///
	keep (1.opower#1.post_opower) ///
	coeflabels(	1.opower 		 		"Opower" ///
				1.opower#1.post_opower 	"Opower $\times$ Post") ///
	label compress nogaps nodepvars ///
	mtitles("All waves" "Wave 3" "Wave 6" "Wave 7") ///
	mgroups("Dependent variable: Electricity Usage in kWh", pattern(1 0 0 0 0 0) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span) ///
	stats(control_mean N_clust N, fmt(%9.0fc) labels("Control mean" "Households" "Observations")) ///
	nonotes addnote("* p<0.10, ** p<0.05, *** p<0.01") ///
	indicate("Demographics = $c" "Baseline usage = $b" `r(indicate_fe)')
	
********************************************************************************
*                   Heterogeneous treatment effect (HTE)                       *
********************************************************************************
use .\intermediate_data\elec_wave367_winterpre.dta, clear

// generate above and below median WITHIN each wave
local covariate "income building_size uni_size member vintage $pre_avg"
foreach var in `covariate' {
	bysort opower_paper_wave: egen `var'_median_wave = median(`var') 
	gen `var'_above_wave = . 
	// arbitrarily use >= for above wave and < for below wave (we do not use binary variables)
	replace `var'_above_wave = 1 if `var'>=`var'_median_wave & `var'!=.
	replace `var'_above_wave = 0 if `var'<`var'_median_wave & `var'!=.
}

rename income_above_wave 		inc_above_wave
rename building_size_above_wave	size_above_wave
rename ${pre_avg}_above_wave 	preuse_above_wave
rename uni_size_above_wave 		unisize_above_wave

// define macros
local above_1 preuse_above_wave
local above_2 size_above_wave
local above_3 inc_above_wave
local above_4 vintage_above_wave
local varlist1 $c
local varlist2 $b $c_b
local varlist3 $b $c_i
local varlist4 $b $c_v


// run regressions
eststo clear
foreach i of numlist 1 2 3 4 {
	cap drop above_wave
	gen above_wave = `above_`i'' // create temporary variable of interest for making consolidated table
	eststo: qui reghdfe tot_kwh above_wave#opower_dt i.opower i.post_opower i.above_wave i.above_wave#opower i.above_wave#post_opower ///
	i.above_wave#c.`varlist`i'' `varlist`i'', absorb(${waveXmonth}#above_wave) vce (cluster ky_ba) 
	test 0.above_wave#1.opower_dt = 1.above_wave#1.opower_dt
	local p : di %9.2f r(p)
	estadd local pval `p'
	qui sum tot_kwh if opower==0 
	local r_mean_hat: di %9.0fc `r(mean)'
	estadd local control_mean `r_mean_hat'
}

* add FE indicators 
estfe . est*, labels(${waveXmonth}#above_wave "Wave $\times$ year-month $\times$ category FE")
return list // create local used in table creation

// Export to tex
esttab using ./output/tables/T4.tex, replace booktabs ///
	b(%9.2f) se(%9.2f) ///
	star(* 0.1 ** 0.05 *** 0.01)  ///
	keep (0.above_wave#1.opower_dt 1.above_wave#1.opower_dt) ///
	coeflabels(	0.above_wave#1.opower_dt "Opower $\times$ Post $\times$ Below Median" ///
				1.above_wave#1.opower_dt "Opower $\times$ Post $\times$ Above Median") ///
	label compress nogaps nodepvars nonumbers ///
	mtitles("Baseline Usage" "House Size" "Income" "House Year Built") ///
	mgroups("Dependent Variable: Electricity Usage in kWh", pattern(1 0) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span) ///
	stats(pval control_mean N_clust N, fmt(%9.0fc) ///
		labels("p-value, test of equal coefficients" "Control mean" "Households" "Observations")) ///
	nonotes addnote("* p<0.10, ** p<0.05, *** p<0.01") ///
	indicate("Demographics $\times$ category = $c_i" `r(indicate_fe)')

********************************************************************************
*                       Summary of estimation sample                           *
********************************************************************************
use .\intermediate_data\elec_wave367_winterpre.dta, clear

bysort opower_paper_wave : gen opower_obs = _N
egen opower_unique = nvals(ky_ba), by(opower_paper_wave)

estpost tabstat opower_obs opower_unique, by(opower_paper_wave) nototal
esttab . using ./output/tables/TA1.tex, replace booktabs ///
	cells("opower_obs(fmt(%9.0fc)) opower_unique(fmt(%9.0fc))") ///
	label noobs nomtitle nonumber ///
	coeflabels(3 "Wave 3" 6 "Wave 6" 7 "Wave 7") ///
	collabels("Number of observations" "Number of unique accounts" )

********************************************************************************
*                           Balance of covariates                              *
********************************************************************************
// conduct balance tests in the cross section
bysort ky_ba (month): gen xs = (_n==1) // generate indicator for cross section

// run ttest and export to tex 
foreach i of numlist 3 6 7 {
	quietly estpost ttest $pre_avg income member building_size uni_size vintage married if opower_paper_wave==`i' & xs==1, by (opower) unequal
	esttab . using ./output/tables/TB_wave`i'.tex, replace booktabs ///
		cells("mu_1(fmt(%09.3gc)) mu_2(fmt(%09.3gc)) b(fmt(%09.3gc) star) t(fmt(2))") ///
		star(* 0.1 ** 0.05 *** 0.01) ///
		label noobs nomtitle nonumber ///
		collabels("Control" "Treatment" "Difference" "t-statistic")
}

********************************************************************************
*                             Robustness check                                 *
********************************************************************************
use .\intermediate_data\elec_wave367_winterpre.dta, clear

eststo clear
clear matrix
local varlist1 $b $c
local varlist2 $b
local varlist3 $c
local varlist4 

foreach i of numlist 1 2 3 4 {
	eststo model0`i': qui reghdfe tot_kwh opower##post_opower `varlist`i'', absorb($waveXmonth) vce(cluster ky_ba)
}

foreach j of numlist 3 6 7{
		foreach i of numlist 1 2 3 4 {
			eststo model`j'`i': qui reghdfe tot_kwh opower##post_opower `varlist`i'' if opower_paper_wave==`j', absorb($waveXmonth) vce(cluster ky_ba)
		}
}

* program written by Ben Jann. We use it to transpose models and append as rows 
program drop appendmodels
program appendmodels, eclass
syntax namelist 
tempname b V tmp 

foreach name of local namelist {
	qui est restore `name'
	mat `tmp' = e(b)
	local eq1: coleq `tmp'
	gettoken eq1 : eq1
	mat `tmp' = `tmp'[1,"`eq1':"]
	local cons = colnumb(`tmp', "_cons")
	if `cons'<. & `cons'>1 {
		mat `tmp' = `tmp'[1,1..`cons'-1]
	}
	mat `b' = nullmat(`b') , `tmp'
	
	mat `tmp' = e(V)
	mat `tmp' = `tmp'["`eq1':","`eq1':"]
	if `cons'<. & `cons'>1 {
		mat `tmp' = `tmp'[1..`cons'-1,1..`cons'-1]
	}
	capt confirm matrix `V'
	if _rc {
		mat `V' = `tmp'
	}
	else {
		mat `V' = ///
		(`V' , J(rowsof(`V'),colsof(`tmp'),0) ) \ ///
		( J(rowsof(`tmp'),colsof(`V'),0) , `tmp')
	}
}

local names: colfullnames `b'
mat coln `V' = `names'
mat rown `V' = `names'
eret post `b' `V' 
eret local cmd "whatever"

end

* compile regression models for transposed pooled 
eststo overall0: appendmodels model01 model02 model03 model04
	qui sum tot_kwh if opower==0  // get pre-treatment mean for pooled sample
	local r_mean_hat: di %9.0fc `r(mean)' // round the number 
	estadd local control_mean `r_mean_hat' // store in local for display in table 
	qui distinct(ky_ba) // get number of observations and number of unique households
	local n_obs = `r(N)'
	estadd local N `n_obs'
	local n_clust = `r(ndistinct)'
	estadd local N_clust `n_clust'

* compile regression models for transposed wave 3, 6, and 7
foreach j of numlist 3 6 7{
eststo overall`j': appendmodels model`j'1 model`j'2 model`j'3 model`j'4
	qui sum tot_kwh if opower==0 & opower_paper_wave==`j'  
	local r_mean_hat: di %9.0fc `r(mean)' 
	estadd local control_mean `r_mean_hat' 
	qui distinct(ky_ba) if opower_paper_wave==`j'
	local n_obs = `r(N)'
	estadd local N `n_obs'
	local n_clust = `r(ndistinct)'
	estadd local N_clust `n_clust'
}
	
esttab overall0 overall3 overall6 overall7 using ./output/tables/T3_appendix_robust.tex, replace booktabs ///
	b(%9.2f) se(%9.2f) ///
	star(* 0.1 ** 0.05 *** 0.01)  ///
	keep(1.opower#1.post_opower) coeflabels(1.opower#1.post_opower 	"Opower $\times$ Post") ///
	refcat(1.opower#1.post_opower "Model specifications in table footnote", nolabel) ///
	mtitles("All waves" "Wave 3" "Wave 6" "Wave 7") ///
	mgroups("Dependent variable: Electricity Usage in kWh", pattern(1 0 0 0 0 0) ///
	prefix(\multicolumn{@span}{c}{) suffix(}) span) ///
	stats(control_mean N_clust N, fmt(%9.0fc) labels("Control mean" "Households" "Observations")) ///
	nonotes addnote("* p<0.10, ** p<0.05, *** p<0.01" ///
	"Row1: including both demographics and baseline usage" ///
	"Row2: excluding demographics" ///
	"Row3: excluding baseline usage" /// 
	"Row4: excluding both demographics and baseline usage") 

********************************************************************************
capture graph close
capture log close
exit
